AD-A234  273 


Computer  Model  of  Atmospheric  Ice 
Accretion  on  Transmission  Lines 


Kathleen  F.  Jones  and  Kurt  Z.  Egelhofer 


February  1991 


•WIND 


100  s 


600  s 


1200  s 


■  '  T  '-5  1 

a  I 


1800  s 


91  4  04  016 


For  conversion  of  SI  metric  units  to  U.S./British  customary  units  of 
measurement  consult  ASTM  Standard  E380,  Metric  Practice 
Guide,  published  by  the  American  Society  for  Testing  and 
Materials,  1916  Race  St.,  Philadelphia,  Pa.  19103. 


Cover:  Cross  section  of  a  transmission  line  with  a 
growing  ice  accretion,  showing  how  the  wire 
rotates  as  the  accretion  gets  larger. 


CRREL  Report  91-3 


U.S.  Army  Corps 
of  Engineers 

Cold  Regions  Research  & 
Engineering  Laboratory 


Computer  Model  of  Atmospheric  Ice 
Accretion  on  Transmission  Lines 

Kathleen  F.  Jones  and  Kurt  Z.  Egelhofer  February  1991 


Prepared  for 

OFFICE  OF  THE  CHIEF  OF  ENGINEERS 


~\ 


\ 


Approved  for  public  release;  distribution  is  unlimited 


PREFACE 


This  report  was  prepared  by  Kathleen  F.  Jones,  Research  Physical  Scientist.  Snow  and  Ice 
Branch,  Research  Division,  U.S.  Army  Cold  Regions  Research  and  Engineering  Laboratory, 
and  Kurt  Z.  Egelhofer,  graduate  student,  Thayer  School  of  Engineering,  Dartmouth  College. 
Funding  for  this  research  was  provided  by  DA  Project  4A 1 6 1 1 02AT24,  Work  Unit  FS/005, 
Concepts  for  Spatial  Winter  Boundary  Layer  Description.  The  authors  thank  Jacqueline 
Richter-Menge  and  Kazuhiko  Itagaki  for  reviewing  the  report. 


CONTENTS 


Page 


Preface .  ii 

Nomenclature .  iv 

Introduction .  1 

Atmospheric  icing  process .  2 

Droplet  trajectories .  2 

Surface  singularity  method .  4 

Finite-element  technique .  6 

Accretion  thermodynamics .  7 

Ice  accretion  profile .  9 

Ice  and  wind  loads .  10 

Computer  model .  12 

Results .  16 

Atmospheric  variables .  17 

Wire  variables .  17 

Summary .  18 

Conclusions .  18 

Future  work .  19 

Literature  cited .  20 

Appendix  A:  Determining  the  angle  of  twist  along  a  cylinder  fixed  at  both  ends  subject 
to  a  distributed  torque  along  its  length .  23 


ILLUSTRATIONS 

Figure 

1 .  Dependence  of  accreted  ice  load  on  the  median  volume  droplet  diameter . 

2.  Air  velocity  field  and  droplet  trajectories . 

3.  Heat  balance  at  the  freezing  surface . 

4.  Horizontal  force  Fx  normalized  to  the  accretion  diameter  DI  vs  the  wind  speed 

normal  to  the  wire  for  rime  ice.  milky  ice  and  clear  ice  accretions . 

5.  Flow  chart  of  the  wire  icing  program . 

6.  Finite-element  mesh  for  the  wire  icing  program . 

7.  Finite-element  mesh  close  to  the  wire . 

8.  Comparison  of  the  local  collection  efficiency  calculated  by  the  present  model 

with  that  of  Lozowski  and  Oleskiw . 

9.  Example  of  wire-ice  surface  temperatures  calculated  by  subroutine  THERMO 

10.  Steps  in  the  calculation  of  the  accreted  ice  layer  profile . 

1 1 .  Final  ice  accretion  profiles  for  1 3  model  simulations  . 

TABLE 

Table 

1.  Results  of  icing  simulations .  16 


6 

7 

11 

12 
13 

13 

14 

15 

16 
17 


iii 


at 


NOMENCLATURE 


A  surface  area  of  heat  transfer  per  meter  of  wire  (m) 

C< j  drag  coefficient 

C’|  lift  coefficient 

Cpa  heat  capacity  of  air  (J/kg  °C) 

Cpc  heat  capacity  per  unit  length  of  wire  and  ice  accretion  (J/m  °C) 
Cpi  heat  capacity  of  ice  (J/kg  °C) 

Cpw  heat  capacity  of  water  (J/kg  °C) 
d  distance  from  center  of  wire  to  ice  center  of  mass  (m) 
dx  horizontal  distance  from  center  of  wire  to  ice  center  of  mass  (m) 
dy  vertical  distance  from  center  of  wire  to  ice  center  of  mass  (m) 

D  diameter  of  wire  or  cylinder  (m) 

D,  diameter  of  wire  and  ice  accretion  normal  to  the  flow  (m) 

£  global  collection  efficiency 

£0  saturation  vapor  pressure  of  water  over  ice  at  T0  (N/m2) 

£s  vapor  pressure  of  water  over  ice  at  Ts  (N/m2) 

£a  fraction  of  water  accreted  as  ice 
£,  force  due  to  weight  of  ice  (N) 

£|  lift  force  (N) 

£v  horizontal  force  (N) 
f  v  vertical  force  (N) 
g  gravitational  acceleration  (m/s2) 

G  shear  modulus  (N/m2) 
h  convective  heat  transfer  coefficient  ( W/m2  °C) 
ice  accretion  thickness  (m) 

./  polar  moment  of  inertia  (m4) 
k  dimensionless  inertia  parameter 
k3  molecular  thermal  conductivity  of  air  (W/m  °C) 

K  wire  torsional  stiffness  (Nm/radian) 

/  wire-ice  surface  vector  at  droplet  impact  point 
L  dimensionless  distance  along  airfoil  surface  from  stagnation  point 
£w  total  length  of  wire  (m) 

£f  latent  heat  of  fusion  of  water  (J/kg) 

£v  latent  heat  of  vaporization  of  water  (J/kg) 
m  ice  accretion  mass  flux  (kg/m2s) 
w,  mass  of  ice  accretion  (kg) 

P  atmospheric  pressure  (N/m2) 

Q  change  in  heat  content  per  meter  of  wire  (J/m) 
r  droplet  radius  (m) 

R  -1  QPrVJTs 
Re  droplet  Reynolds  number 


iv 


Rs  surface  recovery  factor 
t  time  (s) 

r  dimensionless  time 

A t  time  increment  (s) 

T  torque  (Nm) 

Tf  film  temperature  (°C) 

T0  ambient  air  temperature  (°C) 

rs  wire-ice  surface  temperature  (°C) 

U  flow  velocity  (m/s) 

U  dimensionless  flow  velocity 

U0  uniform  far-field  flow  velocity  (m/s) 
ux  x  component  of  flow  velocity  (m/s) 

«v  y  component  of  flow  velocity  (m/s) 

V  droplet  velocity  (m/s) 

l  dimensionless  droplet  velocity 
V0  droplet  impact  speed  at  stagnation  line  (m/s) 

IV  liquid  water  content  (kg/m3) 

v  horizontal  coordinate 

y  vertical  coordinate 

Y  dimensionless  ordinate  of  droplet  trajectory  starting  point 

c  distance  along  wire  (m) 

a  airfoil  angle  of  attack  (radians) 

P  local  collection  efficiency 

y  cylinder  angle,  measured  from  stagnation  line  (degrees) 

8  incremental  rotation  of  wire  (radians) 
stream  function  (m2/s) 

<J>  velocity  potential  (m2/s) 

(j,  dynamic  viscosity  of  air  (kg/m  s) 

v  kinematic  viscosity  of  air  (m2/s) 
pa  density  of  air  (kg/m3) 
p,  density  of  ice  (kg/m3) 
pw  density  of  water  (kg/m3) 

8  rotation  of  wire  (radians) 

0  average  rotation  of  wire  fixed  at  both  ends  (radians) 
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Computer  Model  of  Atmospheric  Ice  Accretion 
on  Transmission  Lines 

KATHLEEN  F.  JONES  AND  KURT  Z.  EGELHOFER 

INTRODUCTION 

Atmospheric  icing  occurs  in  conditions  where  the  cooling  of  an  air  mass  causes  the  supercooling 
of  water  droplets.  Water  droplets  in  the  atmosphere  can  remain  in  the  liquid  state  at  air  temperatures 
as  low  as  -40°C  before  spontaneous  freezing  occurs  (Battan  1962).  Aircraft  operating  in  these 
supercooled  clouds  experience  structural  icing  as  the  water  droplets  hit  the  aircraft  and  freeze.  Icing 
can  occur  at  the  earth's  surface  in  a  supercooled  fog  or  during  freezing  rains.  Under  these  conditions 
ice  will  accrete  on  structures  with  surface  temperatures  below  freezing. 

Ice  accretion  problems  are  of  interest  to  engineers  and  scientists  because  ice  growth  on  manmade 
structures  can  lead  to  structural  failure  and  the  subsequent  loss  of  human  life  or  property.  Widely 
known  examples  of  failure  due  to  ice  growth  include  aircraft  crashes,  tower  collapses  and  transmission 
line  failures  leading  to  power  outages.  Researchers  have  approached  the  accretion  problem  with  three 
techniques:  computer  modeling  studies,  laboratory  studies  and  field  studies.  Recent  attempts  at 
modeling  transmission  line  icing  have  resulted  in  the  quantification  of  many  of  the  important  aspects 
of  the  problem,  including  liquid  water  droplet  interactions  with  the  air  flow,  the  variation  of  local 
collection  efficiency  of  the  accreting  object,  thermodynamic  processes  at  the  freezing  surface,  the 
torsional  response  of  flexible  structures  toasymmetric  loads, and  aerodynamic  forces  on  the  structure. 

Research  on  the  accretion  of  atmospheric  ice  on  ground-based  structures  has  taken  a  variety  of 
directions:  estimating  the  probability  of  an  icing  event,  estimating  the  expected  hazards  of  an  event, 
and  devising  strategies  tomitigate  oreliminate  ice  accretions.  Many  potential  problems  can  be  avoided 
through  prudent  planning  in  conjunction  with  climatological  and  meteorological  studies.  For  example. 
Mallory  and  Leavengood  ( 1983)  described  the  successful  relocation  of  a  Southern  California  Edison 
Company  500-kV  transmission  line  based  on  meteorological  studies  that  predicted  frequent  icing 
events  along  a  segment  of  the  original  route.  However,  atmospheric  icing  is  a  fairly  widespread 
phenomenon,  and  the  expansion  of  activities  into  mountainous  and  subarctic  regions  has  forced 
engineers  and  scientists  to  address  the  problem  directly.  The  lack  of  historical  meteorological  data  in 
remote  locations  has  hindered  attempts  to  determine  icing  parameters  and  to  predict  the  occurrence 
and  severity  of  icing  storms.  For  these  reasons  recent  research  has  concentrated  on  understanding  the 
physical  accretion  process  and  developing  models  to  predict  the  severity  of  icing  events  (for  example. 
Ackley  and  Templeton  1979,  Lozowski  andOleskiw  1983.  McComberetal.  1983,  Smith  and  Barker 
1983).  Model  results  have  been  verified  by  comparison  with  field  and  laboratory  data.  Additionally, 
field  measurements  of  icing  conditions  have  become  more  frequent  in  recent  years  (for  example. 
McComberetal.  1982,  Govoni  and  Ackley  1983,  Krishnasamy  1983,Tattleman  1983).  However,  a 
literature  review  quickly  reveals  how  sparse  the  data  base  is.  making  it  difficult  to  use  existing  data 
for  a  specific  application  or  geographical  area. 

This  study  focuses  on  modeling  the  formation  of  ice  accretions  on  flexible  structures  such  as 
transmission  lines.  The  mode!  requires  the  input  of  meteorological  and  icing  parameters  and  predicts 
the  vertical  and  horizontal  loads  on  the  support  structures  of  the  transmission  line  for  a  given  structural 
design.  The  design  engineer  can  then  decide  to  adjust  the  structural  design  or  the  transmission  line 
location  or  both  to  mitigate  ice  accretion  loads. 


ATMOSPHERIC  ICING  PROCESS 


Ice  accretions  have  been  classified  according  to  source  and  outward  appearance  (Makkonen 
1984a).  Glaze  and  rime  ice,  which  occur  in  many  severe  icing  events,  form  from  supercooled  water 
droplets.  Glaze  is  a  hard,  nearly  bubble-free  and  clear  homogeneous  ice  with  a  density  close  to  0.92 
Mg/m3,  the  density  of  pure  bubble- free  ice.  Glaze  ice  grows  “wet,”  and  water  may  run  back  along  the 
structure.  Rime  ice,  on  the  other  hand,  grows  "dry,”  that  is,  water  droplets  impinging  on  the  structure 
freeze  on  impact.  Hard  rime  is  a  rather  hard,  granular,  white  or  translucent  ice  with  a  density  between 
0.6  and  0.9  Mg/m3.  Soft  rime  is  white  ice  with  a  loose  structure  and  adensity  less  than  0.6  Mg/m3.  There 
are  two  other  types  of  ice  accretions:  hoarfrost,  which  grows  directly  from  water  vapor,  and  wet  snow, 
which  originates  from  snowflakes.  Because  many  severe  icing  events  involve  glaze  and  rime  ice,  this 
study  focused  on  these  types  of  icing. 

Rime  and  glaze  ice  are  formed  when  supercooled  water  droplets  impinge  on  a  structure  and  adhere 
to  the  surface  by  freezing.  Ackley  and  Templeton  (1979)  identified  six  variables  necessary  to  quantify 
the  amount  and  character  of  accreted  ice  on  a  structure:  ambient  air  temperature,  cloud  liquid  water 
content  and  droplet  size  distribution,  wind  speed,  and  cross-sectional  diameter  and  shape  of  the  object. 
Both  the  liquid  water  content  and  the  droplet  size  distribution  are  difficult  to  measure.  Unfortunately 
the  amount  of  ice  that  will  accrete  on  a  wire  is  extremely  sensitive  to  droplet  size,  as  shown  in  Figure 
1  (Makkonen  1984b). The  airtemperature  and  wind  speed  are  more  easily  measured,  but  as  mentioned 
previously,  an  extensive  data  base  for  remote  locations  does  not  exist.  The  shape  of  an  iced  object  may 
be  complex  and  may  change  as  the  ice  accretion  grows  and  the  effective  size  of  the  structure  changes. 
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Figure  1 .  Dependence  of  accreted  ice  load  on  the  median  volume  drop¬ 
let  diameter  (After  Makkonen  1984b.) 


When  a  transmission  line  is  exposed  to  icing  conditions,  the  rate  of  ice  accretion  is  governed  by  two 
processes  that  depend  on  the  above  variables:  the  impingement  of  supercooled  water  droplets  on  the 
wire,  and  the  thermodynamics  at  the  wire  surface,  which  determines  what  portion  of  the  impinging 
water  freezes  or,  on  the  other  hand,  melts  previously  accreted  ice.  Each  of  these  processes  v  ill  be 
discussed  from  the  standpoint  of  how  previous  researchers  have  approached  the  problem. 


DROPLET  TRAJECTORIES 

For  flow  velocities  at  which  atmospheric  icing  occurs,  the  air  flow  around  a  long  uniform  object 
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can  be  treated  as  two-dimensional,  steady,  incompressible,  irrotational  fluid  flow.  Therefore,  potential 
theory  can  be  used  to  describe  the  velocity  field,  and  the  flow  satisfies  Laplace’s  equation  in  the  entire 
plane  (Batchelor  1970): 

d-(J>/d.i-  +  cP<! >/3>~  =  0  +  3-xjr/dy2  =  0  ( 1 ) 


where.v  =  horizontal  coordinate 
y  =  vertical  coordinate 
0  =  velocity  potential 
\y  =  stream  function. 

Lines  of  constant  y  are  the  streamlines  of  the  flow,  which  are  perpendicular  to  the  lines  of  constant 
(J).  The  air  velocity  components  can  be  determined  at  any  point  in  the  air  stream  by  differentiating  the 
stream  function: 

ux  =  d\\l/dy  Uy  =  -d\\i/d x  (2) 

where  is  the  ,v  component  of  the  flow  velocity  and  uvis  the  y  component  of  the  flow  velocity.  The 
stream  function  and  velocity  potential  are  related  through  the  Cauchy-Riemann  equations  of  complex 
variable  theory: 

dy/3y  =  3<f)/3.v  3y/ck  =  -d0/3v.  (3) 

For  cylinders  with  circular  or  elliptical  cross  sections,  the  velocity  at  any  point  in  the  airstream  can 
be  determined  analytically  (Ackley  and  Templeton  1979).  For  arbitrary  cross-sectional  shapes,  such 
as  airfoils  or  iced  transmission  lines  that  are  free  to  rotate,  transformation  or  numerical  techniques  must 
be  employed.  These  techniques  can  be  broadly  classified  into  three  groups: 

•  Conformal  transformation  techniques  (Theodorson  and  Garrick  1932); 

•  Surface  singularity  methods  (Fless  and  Smith  1967);  and 

•  Finite-element  techniques  (McComberandTouzot  1981). 

Lozowski  and  Oleskiw  (1983)  emphasized  that  the  technique  chosen  should  be  efficient  in  terms  of 
computer  time  because  of  the  large  number  of  velocity  calculations  necessary  to  compute  the  water 
droplet  trajectories,  and  it  should  be  capable  of  adjusting  to  a  changing  object  cross  section  as  ice 
accretes. 

To  determine  the  kinematic  interaction  between  the  water  droplets  in  the  air  stream  and  a  trans¬ 
mission  line,  the  equation  of  motion  for  the  droplets  must  be  solved.  Ackley  and  Templeton  (1979) 
based  their  analytic  approach  for  a  cyl  inder  with  an  elliptical  cross  section  on  the  work  of  Brun  (1957). 
T.  assumptions  for  this  formulation  are: 

•  The  concentration  of  droplets  in  the  air  stream  is  sufficiently  small  that  the  flow  is  not  perturbed 
by  the  droplets; 

•  The  force  of  gravity  is  much  less  than  the  inertial  forces  and  can  be  neglected;  and 

•  The  pressure  forces  on  a  droplet  are  equal  to  those  on  an  equivalent  volume  of  air  at  the  same 
location,  and  these  forces  can  be  neglected  because  the  density  of  water  is  much  greater  than  that 
of  air. 

The  motion  of  the  droplet  therefore  depends  primarily  on  its  inertia  and  the  viscous  drag  force  on  the 
droplet  due  to  the  deflection  of  the  air  stream  around  the  cylinder. 

Langmuir  and  Blodgett  ( 1 946)  formulated  this  problem  for  a  cylinder  and  solved  Newton's  second 
law  of  motion  in  dimensionless  form: 

k  dV/dr  =  -Cd Re{V-u)/ 24  (4) 
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where  Cd  = 
V  = 
U  = 

t  = 

Re  = 
k  = 
= 

(7  = 
V'  = 

I  = 

Pw  = 
Pa  = 
P  = 
r  = 
D  = 


drag  coefficient 

v/u0 

U!U0 

-tUJD 

2pa/-|  V-U\  /|i  (droplet  Reynolds  number) 

Ap^r-U J(A\xD)  (inertia  parameter) 

uniform  far-field  flow  velocity 

flow  velocity 

droplet  velocity 

time 

density  of  water 
density  of  air 
dynamic  viscosity  of  air 
droplet  radius 
cylinder  diameter. 


The  functional  dependence  of  the  drag  coefficient  of  a  spherical  droplet  in  an  air  stream  on  the 
Reynolds  number  has  been  determined  experimentally  by  several  researchers.  Beard  and  Pruppacher 
(1969)  leported  the  following  formulas  for  Cd: 

Cd  =  24  ( 1  +  0. 1 02  Re0  955)//?e  for  0.2  <  Re  <  2 


=  24(1  +0.115  Re0M2)/Re 


for  2  <  Re  <  21 


(5) 


=  24(1  +0.189  Reil6i2)/Re 


for  21  <  Re  <200. 


Previous  researchers,  including  Ackley  and  Templeton  (1979),  solved  eq  4  by  numerical  inte¬ 
gration.  They  then  calculated  a  global  collection  efficiency  E,  which  is  the  fraction  of  water  droplets 
in  the  path  of  the  cylinder  that  collide  with  it.  The  global  collection  efficiency  is  calculated  by  finding 
the  droplet  trajectory  that  is  tangent  to  the  cylinder  and  determining  the  initial  far-field  distance  from 
the  centerline  of  the  cylinder  axis  of  that  tangent  trajectory.  For  example,  if  the  initial  far-field  ordinate 
of  the  tangent  trajectory  is  0.75  units  (with  the  cylinder  radius  representing  one  unit),  then  E  =  0.75. 
Foradistribution  of  waterdroplet  sizes,  the  total  collection  efficiency  is  calculated  by  taking  the  global 
collection  efficiency  for  each  droplet  size,  multiplying  by  the  fraction  of  the  total  liquid  water  content 
represented  by  that  size,  and  summing  over  all  droplet  sizes  in  the  distribution. 

Solving  the  equation  of  motion  in  this  manner  allows  time-dependent  effects  associated  with  the 
droplet  trajectories  to  be  included  in  the  model.  As  ice  accretes,  the  object  changes  shape  and  the  flow 
field  around  the  object  changes.  Ackley  and  Templeton  ( 1 979)  modeled  this  effect  by  conforming  the 
accretion  to  an  elliptical  shape.  Since  an  analytical  solution  to  the  potential  flow  around  an  elliptical 
cylinder  exists,  they  could  determine  the  air  flow  field  analytically.  The  drawback  to  this  formulation, 
however,  is  that  it  is  applicable  only  to  objects  with  this  simple  cross  section.  Using  numerical  and 
transformation  techniques,  several  researchers  have  developed  models  that  can  be  applied  to  objects 
of  arbitrary  cross  section,  including  airfoils  or  wires  that  rotate  freely.  Two  possible  approaches  are 
the  surface  singularity  method  and  the  finite-element  technique. 

Surface  singularity  method 

Lozowski  and  Oleskiw  (1983)  developed  a  computer  model  for  airfoil  icing  that  permitted 
simulation  of  the  time-dependent  growth  of  ice  without  runback  on  an  arbitrary  two-dimensional 
airfoil.  The  model  can  be  used  to  predict  the  icing  rate  as  well  as  accretion  shape  on  airfoils  when  the 
accretion  is  dry.  Dry  accretion  occurs  when  all  of  the  impinging  water  droplets  freeze  upon  impact  so 
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there  is  no  runback  icing.  This  requires  that  the  heat  transfer  from  the  airfoil  surface  is  sufficient  to 
freeze  all  of  the  impinging  water  droplets  while  maintaining  the  surface  temperature  below  03C. 

Lozowski  and  Oleskiw  used  a  surface  singularity  method  to  determine  the  flow  field  around  the 
two-dimensional  airfoil.  This  technique  divides  the  airfoil  surface  into  a  series  of  straight-line 
segments,  and  a  constant  but  unknown  vorticity  is  distributed  along  each.  The  stream  function  at  any 
point  external  to  the  airfoil  is  then  calculated  by  solving  a  set  of  linear  algebraic  equations  for  the 
unknown  vorticity  densities.  The  air  velocity  components  are  calculated  at  any  desired  point  in  the  air 
stream  by  numerically  differentiating  the  stream  function.  This  approach  has  an  advantage  overother 
techniques  in  that  the  velocity  components  are  determined  at  every  point  along  the  droplet  trajectory, 
rather  than  being  interpolated  from  a  velocity  field  calculated  at  a  discrete  array  of  grid  points. 
Lozowski  and  Oleskiw  then  applied  the  droplet  equation  of  motion,  which  describes  the  muiion  of  a 
spherical  droplet  in  an  accelerated  air  flow: 

dV  _  2(Pw~Pa),g  3CdpjV-^(V-t/)  9paV(v/it)  p  dt/tft  (&) 

*  (2Pv.+Pj  4'(2PW+Pa)  (2pw+Pa)  L  drV(t-T) 

buoyancy  drag  history 

where  g  is  gravitational  acceleration  and  v  is  the  kinematic  viscosity  of  air.  The  buoyancy  term 
describes  the  vertical  acceleration  of  the  droplets  due  to  gravity.  The  history  term  takes  into  account 
the  fact  that  droplets  are  accelerating  with  respect  tothe  flow  and  that,  in  fact,  the  problem  is  not  steady- 
state.  This  equation  also  incorporates  the  induced  droplet  mass  that  results  from  the  acceleration  of  air 
in  the  vicinity  of  an  accelerating  droplet.  The  buoyancy  and  history  terms  and  induced  mass  effects 
are  frequently  dropped  in  many  ice  accretion  models,  resulting  in  the  simplified  form  of  this  equation 
that  was  presented  as  eq  4.  Lozowski  and  Oleskiw  soivedeq  6  numerically  using  a  fourth-order  Runge- 
Kutta-Fehlberg  method,  with  the  history  term  approximated  by  a  combined  numerical  and  analytical 
technique.  Although  their  tests  were  not  exhaustive,  they  indicated  that  the  history  tend  could  be 
omitted  without  signified. riy  affecting  the  results.  The  history  term  is  most  significant  in  cases  with 
low  collection  efficiencies. 

Lozowski  and  Oleskiw  used  local  collection  efficiencies  todetermine  the  rate  of  ice  growth  on  each 
small  segment  of  the  airfoil  surface.  The  advantage  of  this  approach  over  earlier  models  (e.g.  Ackley 
and  Templeton  1979).  which  use  a  global  collection  efficiency  to  calculate  the  rate  of  ice  growth  in 
a  predetermined  shape,  is  its  ability  to  model  ice  growth  on  objects  with  an  arbitrary  shape.  The  local 
collection  efficiency  P(Z.)  at  any  point  on  the  airfoil  surface  is  given  by 

P(L)  =  dY/dL  (7) 

where  L  is  the  distance  along  the  airfoil  surface  from  the  stagnation  point  and  Y  is  the  ordinate  of  the 
droplet  starting  point.  The  variables  used  in  eq  7  are  depicted  in  Figure  2.  Lozowski  and  Oleskiw 
calculated  the  trajectories  of  10-20  droplets  with  different  starting  ordinates  and  fit  the  resulting  Y.L 
pairs  to  a  cubic  spline,  which  is  differentiated  toobtain  the  local  collection  efficiency.  They  calculated 
the  accretion  thickness  /r,  by 

h,  =  m{L)M  pj  (g) 

where  iii(L)  =  U0W$(L)  (9) 

Af  =  time  period  of  accretion 
pi  =  ice  density 
W  =  liquid  water  content. 

Using  these  equations  a  new  airfoil  surface  shape  is  determined  for  each  specified  time  interval. 
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Figure  2.  Air  velocity  field  (arrows)  and  droplet  trajectories  (solid 
lines).  Also  shown  are  dL  and  d Y  used  in  eq  7.  (After  Lozowski  and 
Oleskiw  1 9S3 . 1 


Finite-element  technique 

Finite-element  models  employ  a  piecewise  approximation  of  the  solution  to  the  governing 
differential  equations  over  the  domain  of  interest.  The  primary  advantage  of  this  method  is  its  ability 
to  model  irregularly  shaped  boundaries,  even  when  the  boundary  is  changing.  Another  advantage  is 
that  the  size  and  shape  of  the  finite  elements  can  be  varied,  allow  ing  greater  accuracy  in  specific  areas 
of  the  domain  by  using  more  elements  of  smaller  size  there.  The  finite-element  method  does  have  some 
dixads  antages.  The  finite-element  grid  must  be  designed  with  the  expected  solution  to  the  problem  in 
mind.  For  example,  in  this  problem,  in  regions  where  the  air  flow'  velocity  gradients  are  large,  the 
elements  must  be  relatively  small.  Also,  since  the  grid  can  be  irregular,  establishing  the  coordinates 
of  each  node  can  be  tedious.  This  difficulty  can  be  overcome  by  using  a  mesh-generating  routine. 

The  basic  steps  that  must  be  taken  when  using  the  finite-element  method  are  as  follows.  First,  the 
region  of  interest,  in  this  case  the  air  in  front  of  the  wire,  must  bedivided  into  elements,  either  triangles 
or  quadrilaterals.  Next,  the  interpolating  function  to  be  used  it  eh  element  is  chosen.  Polynomials 
ttre  generally  used  because  they  can  be  easily  differentiate^  and  integrated.  The  third  step  is  the 
derivation  of  the  equations  that  describe  the  properties  of  each  element.  These  equations  are  developed 
in  matrix  form,  w  ith  the  size  of  the  matrix  equal  to  the  number  of  nodes  in  each  element.  In  many 
simulations,  as  w  ith  our  model,  Galerkin's  formulation  is  used.  The  fourth  step  is  to  assemble  the 
element  equations  into  a  single  global  matrix  of  equations  that  describes  the  entire  region  of  interest. 
This  global  matrix  is  then  modified  to  account  for  the  boundary  conditions  of  the  problem.  The  result 
is  a  system  of  simultaneous  equations,  with  the  number  of  equations  equal  to  the  number  of  nodes  in 
the  problem.  Finally,  the  system  of  equations  is  solved  using  standard  matrix  decomposition  pro¬ 
cedures. 

McComber  (1981)  and  McContber  and  Touzot  (1981)  developed  a  numerical  model  of  ice 
accretion  using  the  finite-element  method.  The  air  and  droplet  velocity  fields  were  calculated  using 
triangular  isoparametric  quadratic  elements.  The  finite-element  method  was  chosen  because  the  gnd 
could  be  altered  easily  for  different  shapes.  The  droplet  trajectories  w'ere  determined  by  solving  eq  4. 
The  drag  coefficient  of  the  impinging  water  droplets  was  calculated  using  Beard  and  Pruppacher's 
( 1969)  formulation,  shown  previously  as  eq  ?.  McComber  ( 1981 )  then  u..-'dthe  droplet  velocity  field 
to  determine  the  local  collection  efficiency  (3  on  the  cylinder  from  the  following  relationship: 

P=  I'  d// 1  d/I  ,10) 


where  V  is  the  dimensionless  droplet  velocity  at  impact  and /is  the  surface  direction  vector  althe  impact 
point.  The  mass  flux  of  droplets  impinging  on  the  surface  of  the  cylinder  was  determined  by  taking 
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the  product  of  the  local  collection  efficiency,  the  liquid  water  content,  and  the  free  stream  velocity. 
Their  simulation  was  for  dry  growth  conditions  only,  and  thus  a  thermodynamic  balance  at  the  surface 
of  the  cylinder  was  not  required.  Model  results  were  compared  with  experimental  results  from  wind 
tunnel  tests,  and  qualitative  agreement  was  found.  However,  the  numerical  model  consistently 
overestimated  the  amount  of  ice  accreted  on  the  cylinders.  Lozowski  et  al.  ( 1 985)  commented  on  the 
definition  of  the  collection  efficiency  shown  in  eq  10  and  indicated  that  it  is  incorrect  and  tends  to  give 
larger  collection  efficiencies  than  actually  occur.  This  may  explain  the  overprediction  of  ice  accretion 
mass  in  McComber  (1981). 


ACCRETION  THERMODYNAMICS 


Once  the  collection  efficiency  of  the  wire-ice  system  has  been  found,  the  next  task  is  to  determine 
whether  the  droplets  freeze  on  impact.  Several  researchers  (McComber  1981,  Lozowski  and  Oleskiw 
1 983)  assumed  dry  growth  conditions  in  their  models,  in  which  all  of  the  impinging  water  freezes  at 
the  point  of  impact.  This  is  a  good  assumption  at  lower  liquid  water  contents  and  lower  air  temperatures 
when  the  latent  heat  released  during  freezing  is  removed  from  the  wire  more  quickly  than  it  is  released 
by  the  droplets  freezing.  However,  at  higher  liquid  water  contents  and  at  temperatures  approaching 
0°C,  wet  growth  conditions  frequently  exist. 


Aerodynamic  Heating  (  +  ) 


Figure  3.  Heat  balance  at  the  freezing  surface.  (From  Acklev  and  Templeton 
1979.) 


Modeling  wet  accretions  is  more  difficult  because  only  a  fraction  of  the  impinging  water  freezes 
on  impact.  The  remaining  water  runs  back  toward  the  separation  region  of  the  wire,  where  it  freezes 
oris  shed.  Ackley  and  Templeton  ( 1 979)  calculated  the  heat  balance  at  the  surface  of  a  cylinder  during 
ice  growth,  basing  their  calculation  on  three  conditions  defined  by  whether  all,  a  fraction,  or  none  of 
the  water  froze  on  impact.  Figure  3  shows  the  heat  balance,  indicating  the  direction  of  heat  flow.  The 
heat  balance  is  the  sum  of  the  aerodynamic  heating  from  compression  of  the  flow,  the  kinetic  energy 
of  droplet  impact,  the  latent  heat  released  by  droplet  freezing,  the  convective  and  evaporative  heat  lost 
to  the  air  flow,  and  the  heat  lost  in  warming  supercooled  droplets  to  0°C: 


I  2 

Q  =  -hAAt !  (7V-7V)  -  + 

*  ~>c 

L  pa 

convective  aerodynamic 


0.6Ev(Es-Eo) 

evaporative 


-AmCpwAf 


1  Cpw 


T  \  Vt  Up 

1  Cpw  2Cpw 


supercooling  latent  kinetic 


(11) 
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where  Q  = 
h  = 
A  = 


r  = 

'~pl.w 

Li  = 


A/  = 
= 


change  in  heat  content  per  meter  of  wire 

heat  transfer  coefficient 

surface  area  of  heat  transfer  per  meter  of  wire 

wire-ice  surface  temperature 

ambient  temperature 

surface  recovery  factor  (0  <  Rs  <  1 ) 

heat  capacity  of  air 

latent  heat  of  vaporization  of  water 

saturation  vapor  pressure  of  water  over  ice  at  Ts 

saturation  vapor  pressure  of  water  over  ice  at  Ta 

atmospheric  pressure 

heat  capacity  of  water 

heat  capacity  of  ice 

Cpi  for  Ts  <  0  or  Cpw  for  Ts  >  0 

latent  heat  of  fusion  of  water 

time  increment 

fraction  of  water  accreted  as  ice. 


In  dry  growth  conditions,  Ts  <  0°C,  Fa  =1  and  Cpj  is  used  in  the  supercooling  term.  In  wet  growth,  7S 
=  0°C  and  Fa  is  between  0  and  1 .  For  77  >  0°C  none  of  the  water  accretes  as  ice,  so  Fa  =  0  and  Cpw  is 
used  in  the  supercooling  term.  The  heat  transfercoefficient  is  taken  from  Bosch's  formulation  (Boelter 
etal.  19651: 

/t  =  31.01  ",f> [273. \hP  /  (101 4T,)ja?6 /D0'44  ( 1 2) 

where  7j  =  0.5(7o-tTs). 

In  their  analysis  Ackley  and  Templeton  (1979)  followed  Messinger  (1953)  and  applied  the 
appropriate  temperature  condition  to  solve  for  the  change  in  heat  content  of  the  cylinder-ice  system, 
which  is  then  used  to  adjust  the  surface  temperature  and  heat  capacity  Cpc  of  the  wire  and  accretion 
during  a  time  interval  Ar  using  the  following  equations: 


K  =  T^Q/Cpc 

(13) 

C’pc=Cpc+FfPimAL<  • 

(14) 

Terms  with  a  prime  denote  updated  quantities.  For  Ts  =  0°C  the  fraction  accreted  is  initially  assumed 
to  be  unity,  and  if  the  resulting  heat  transfer  is  sufficient  to  freeze  all  of  the  water,  the  temperature  drops 
below  0°C  and  dry  accretion  conditions  apply.  If  not,  an  iterative  procedure  is  initiated  to  determine 
the  fraction  accreted  for  the  calculated  change  in  heat  content. 

Other  factors  could  be  included  in  the  thermodynamic  model.  Giedt  (1949)  showed  from  ex¬ 
perimental  work  that  the  convective  heat  transfer  varies  by  a  factor  of  about  two  between  the  stagnation 
line  and  separation  points  on  acircular  cylinder  fora  Reynold's  number  based  on  thecylinderdiameter 
ranging  from  7  x  104  to  2.2  x  I05.  Lozowski  et  al.  ( 1979)  developed  a  model  for  icing  with  runback 
of  unfrozen  water  on  a  fixed  cylinder.  In  their  model  they  used  a  convective  beat  transfer  coefficient 
that  varied  with  angle  y  measured  around  the  cylinder  from  the  stagnation  line.  For  a  smooth  cylinder 

h  =  *a|t/„ Pa/(2n)|0-5  [  1  -  (y/90o)’l  (15) 


where  k,d  is  the  molecular  thermal  conductivity  of  air.  This  relationship  is  based  on  the  data  of 
Achenbach  ( 1974)  and  is  valid  for  y  between  0°  and  about  80°.  Another  variable  that  affects  heat 
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transfer  is  the  roughness  of  the  accretion  surface.  Preliminary  measurements  by  Achenbach  (1977) 
indicated  that  the  convective  heat  transfer  coefficient  on  a  rough  cylinder  may  be  as  much  as  twice  that 
of  a  smooth  one. 


ICE  ACCRETION  PROFILE 

The  size  and  shape  of  the  ice  accretion  are  determined  by  the  density  of  the  accreting  ice  and  the 
torsional  flexibility  of  the  wire.  The  eccentric  weight  of  the  accreted  ice  will  cause  the  flexible  wire 
to  rotate  as  the  ice  accretes.  According  to  Macklin  ( 1 962)  the  average  density  p,  of  rime  ice  accreted 
on  a  rotating  cylinder  depends  on  the  droplet  size,  the  droplet  impact  velocity  at  the  stagnation  point, 
and  the  surface  temperature  of  the  accretion  and  is  given  by  the  following  relationship: 

Pi  =  0.1 10/?°  7(1  Mg/m3  (16) 

where  R  -  -1 06rV'Q/7's  and  T0  is  the  droplet  impact  speed  at  the  stagnation  line.  This  equation  applies 
for/?  values  between  0.8  and  10;  Macklin  suggested  taking  Pj  =  0.9  Mg/m3  for  larger  values  of  R.  Bain 
and  Gayet(1983)  provided  a  formula  that  approximates  Macklin's  curve  for/?  between  10  and  60; 

pi  =  /?/(/?  + 5.61 )  Mg/m3.  (17) 

For  larger  values  of  /?,  they  took  the  accreted  ice  density  to  be  0.9 1 7  Mg/m3. 

As  ice  accretes  on  the  windward  side  of  a  transmission  line,  the  wire  rotates  in  response  to  the 
eccentric  weight.  The  average  incremental  rotation  Sofa  wire  due  to  a  total  applied  torque  T  is 

T  =  K(6+5)  (18) 


where  K  is  the  average  wire  torsional  stiffness  and  9  is  the  wire  rotation  due  to  previous  torque  on  the 
wire.  Following  Smith  and  Barker  (1983)  the  torque  on  a  wire  due  to  an  eccentric  mass  mt  is 

T  =  m\gcl  sin[tan-l(</x/</y)  +  5]  ( 19) 


w  here  dx  =  horizontal  distance  from  the  center  of  mass  of  the  ice  accretion  to  the  wire  center 
<4,  =  vertical  distance  from  the  center  of  mass  of  the  ice  accretion  to  the  wire  center 


At  equilibrium  these  torques  will  be  equal.  Assuming  5  is  small  (say.  less  than  0.1  radian),  then 


5  =  (/np'dj-X'B)  /  iK-ni^’cly).  (20) 

Thus,  the  incremental  rotation  of  a  wire  depends  on  the  stiffness  and  previous  rotation  of  the  wire  as 
well  as  the  mass  and  moment  arm  of  the  accretion. 

The  effective  torsional  stiffness  of  a  wire  fixed  at  both  ends  varies  with  distance  from  the  fixed  ends. 
In  our  model  the  wire  rotation  is  calculated  using  the  average  torsional  stiffness  over  the  entire  length 
of  the  wire.  The  rotation  of  a  wire  fixed  at  both  ends  and  subject  to  a  uniform  torque  T  along  its  length 
Fw  is 


9(r)  =  TLw(z/Lw)  [  M.-/LW))/(2G7) 


(21) 


where  G  =  wire  shear  modulus 

J  -  wire  polar  moment  of  inertia 


9 


z  =  distance  along  the  wire 
Lw  =  wire  length. 


This  equation  is  derived  in  Appendix  A.  The  average  stiffness  can  be  found  by  determining  the  average 
rotation  over  the  length  of  the  wire: 

0  =  -L  [%  (;  )dz  =  TL  w/(  1 2GJ ).  (22) 

Lw  0 

This  gives  an  average  wire  stiffness  of 

K=\2GJILV..  (23) 

According  to  McComber  ( 1983),  for  a  typical  aluminum  stranded  cable  (Bersimis  conductor  with  D 
=  3.5  cm),  GJ  =  351.4  Nm2/rad.  A  250-m-long  wire,  for  example,  will  have  an  average  torsional 
stiffness  of  16.9  Nm/rad. 

ICE  AND  WIND  LOADS 

The  ultimate  reason  for  modeling  the  ice  growth  on  a  transmission  line  during  an  icing  event  is  to 
predict  the  horizontal  and  vertical  forces  exerted  on  the  wire  and  its  supporting  structures.  These  forces 
are  due  to  the  combined  effects  of  the  weight  of  ice  on  the  wire  and  any  wind  loading  either  during  or 
following  an  icing  event.  Transmission  line  support  structures  are  generally  designed  fora  vertical 
compressive  load  imposed  by  the  weight  of  the  conductor  and  any  accreted  ice.  The  horizontal  force 
due  to  the  wind  drag  on  an  iced  wire  is  difficult  to  quantify  because  it  is  a  function  of  the  size,  shape 
and  roughness  of  the  accretion.  Currently  available  models  are  capable  of  predicting  the  size  and,  in 
some  cases,  the  shape  of  the  accretion  but  not  the  roughness,  which  depends  on  the  growth  regime.  The 
numerical  simulation  by  Makkonen  (1984b)  indicates  that  the  growth  regime  may  change  from  wet 
growth  to  dry  growth  during  an  icing  event  of  constant  atmospheric  conditions. 

Assuming  pseudostatic  conditions,  the  total  horizontal  force  due  to  wind  drag  on  a  length  of  wire 
with  a  circular  cross  section  is  (Blevins  1984) 

(24) 

where  Cd  =  drag  coefficient 
pa  =  density  of  air 

Dj  =  diameter  of  the  wire  and  ice  accretion  perpendicular  to  the  airflow  direction. 

For  smooth  cylinders  the  drag  coefficient  is  1 .2  at  low  Reynolds  numbers  and  drops  as  low  as  0.4  as 
the  Reynolds  number  increases  and  the  boundary  layer  becomes  turbulent.  However,  ice  accretions 
are  typically  neither  cylindrical  nor  smooth.  Assuming  that  the  effects  of  roughness  and  shape  may 
offset  each  other  somewhat  (the  overall  shape  may  be  somewhat  streamlined,  and  the  increase  in  drag 
due  to  roughness  may  be  partially  offset  by  the  drop  in  drag  coefficient  associated  with  the  transition 
to  a  turbulent  boundary  layer),  then  taking  Cd=  1  may  not  be  unreasonable.  Similarly  the  vertical  force 
due  to  lift  on  the  wire  is  (Blevins  1984) 

F,  =  0.5Cfi^yo  (25) 

where  Ct  is  the  lift  coefficient.  For  airfoil  cross  sections  the  lift  coefficient  is  proportional  to  the  angle 
of  attack  a,  for  a  less  than  approximately  1 5°.  As  the  absolute  angle  of  attack  increases,  the  point  of 
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Figure  4.  Horizontal  force  Fv  normalized  to  the  accretion  diameter 
D,  vs  the  wind  speed  normal  to  the  wire  for  rime  ice,  milky  ice  and 
clear  ice  accretions.  The  hest-fit  line  for  the  milky  ice  accretion  is 
shown.  (After  Govoni  and  Ackley  1983.) 

flow  separation  moves  forward  until  a  critical  angle  of  attack  is  reached  w  here  stall  occurs.  Then  the 
magnitude  of  the  lift  coefficient  decreases  rapidly  and  the  drag  force  increases. 

The  U.S.  Army  Cold  Regions  Research  and  Engineering  Laboratory,  in  conjunction  w  ith  the  Mt. 
Washington  Observatory,  has  been  collecting  simulated  transmission  line  icing  data  since  1 978  at  the 
top  of  Mt.  Washington,  New  Hampshire  (Govoni  and  Ackley  1983).  Instrumentation  for  the  study 
consisted  of  a  triaxial  load  cell  mounted  at  one  end  of  an  8-m  wire  capable  of  measuring  vertical, 
horizontal  and  axial  forces.  They  observed  that  the  wire  twisted  as  rime  formed  on  its  windward  side, 
a  consequence  of  the  asymmetric  weight  of  the  accretion.  As  a  result  the  final  cross-sectional  structure 
of  the  ice  formation  had  a  spiral  pattern. 

Govoni  and  Ackley  (1983)  analyzed  the  relationship  between  the  drag  force  on  the  ice  accretion, 
normalized  to  the  diameter  of  the  accretion,  and  the  wind  speed.  Figure  4  shows  a  plot  of  FJDt  vs  U0 
for  milky  ice  (high  density  with  intermediate  roughness),  clear  ice  and  rime.  The  plot  indicates  a  good 
correlation  between  the  normalized  drag  force  and  the  wind  speed  for  milky  ice,  implying  a  relatively 
constant  drag  coefficient  throughout  the  icing  event.  However,  when  all  the  ice  types  are  included,  the 
correlation  between  the  normalized  drag  force  and  the  wind  speed  is  weak,  which  was  attributed  to 
changes  in  the  drag  coefficient  due  to  the  changing  shape  of  the  accretion  as  ice  collected.  Figure  4 
does  indicate,  however,  an  increase  in  the  normalized  drag  force  at  a  given  wind  speed  for  increasingly 
rough  accretions. 

According  to  McComberet  al.  ( 1 983),  Hydro-Quebec  transmission  lines  are  designed  to  withstand 
three  loads:  a  wind  load  on  bare  wires  corresponding  to  a  wind  speed  of  1 37  km/hr  (85  mph);  an  ice 
load  consisting  of  a  maximum  radial  ice  thickness  of  4.45  cm  ( 1 .75  in.)  with  no  wind:  and  a  combined 
load  corresponding  to  a  wind  speed  of  72  km/hr  (45  mph)  and  a  radial  ice  thickness  of  3. 18  cm  ( 1 .25 
in.).  The  following  method  is  used  to  calculate  these  design  loads.  The  horizontal  force  due  to  wind 
drag  on  a  bare  wire  is  determined  from  eq  24,  with  Cd  =  I  and  pa  =  1.3  kg/m3  to  give 

Fx  =  0.65  DL^ul  (26) 

The  horizontal  force  due  to  the  combined  effects  of  wind  and  ice  are  calculated  using  eq  26.  with  the 
exception  that  the  ice  accretion  diameter  D,  is  used  for  D.  This  diameter  is  calculated  from  the  ice 


weight,  assuming  an  ice  density  of  0.9  Mg/m\  The  vertical  force  Fv  is  obtained  by  adding  the  force 
due  to  the  weight  of  accreted  ice  Ft  and  the  force  due  to  the  weight  of  the  wire.  McComber  et  al.  ( 1 983) 
conducted  wind  tunnel  experiments  in  which  soft  rime,  hard  rime  and  glaze  were  formed  on  rods. 
Horizontal  and  vertical  forces  on  the  iced  rods  were  measured  as  a  function  of  wind  speed.  In  all  cases 
the  measured  forces  were  about  twice  as  big  as  the  forces  they  calculated.  It  was  determined  that  the 
measured  horizontal  forces  were  greater  for  two  reasons:  a)  the  actual  rough  asymmetric  shapes  of  the 
ice  accretions  had  drag  coefficients  larger  than  the  assumed  Cd  =  1  and  b)  the  soft  and  hard  rime 
accretions  had  lower  densities  than  the  assumed  0.9  Mg/nv'  and  therefore  had  larger  diameters  than 
were  calculated  from  their  weights,  resulting  in  a  further  underestimate  of  the  drag  force.  Also,  it  was 
found  that  in  all  cases  except  one,  the  mea¬ 


sured  lift  force  was  negative  and  thus  was 
equivalent  to  additional  weight.  They  con¬ 
cluded  that  it  will  be  difficult  to  predict 
aerodynamic  forces  on  a  power  line  because 
the  final  asymmetric  shape  of  the  accretion 
depends  strongly  on  the  history  of  its  forma¬ 
tion.  The  results  of  this  study  support  the 
development  of  an  ice  accretion  model  to 
predict  icing  loads  on  wires  with  finite  tor¬ 
sional  rigidities. 

COMPUTER  MODEL 

This  study  draws  from  previous  icing 
models  toconstruct  an  icing  model  fora  wire 
of  an  arbitrary  non-zero  stiffness,  allow  ing 
the  size  and  shape  of  the  ice  accretion  to  be 
determined  locally  by  the  collection  effi¬ 
ciency.  the  surface  heat  balance,  the  ice 
density  and  the  wire  stiffness.  The  model 
consists  of  a  main  program,  w  here  input  and 
output  are  manipulated  and  all  decisions 
occur,  and  thirteen  subroutines,  where  all 
calculations  are  made.  Figure  3  is  a  flow 
chart  of  the  model,  with  each  box  represent¬ 
ing  a  subroutine.  In  the  first  part  of  the 
model,  the  water  droplet  trajectories  and  the 
wire  collection  efficiency  for  a  particular 
droplet  size  are  determined.  In  the  second 
part  of  the  model,  the  heat  balance  calcula¬ 
tions  are  carried  out.  the  accretion  profile  is 
determined,  and  the  rotation  of  the  wire  due 
to  the  asymmetric  ice  accretion  iscalculated. 

The  model  starts  by  reading  the  model 
input,  which  consists  of  variables  that  define 
the  problem  to  be  solved,  including  droplet 
size,  ambient  air  temperature,  liquid  water 


content,  wind  speed,  wire  diameter  and  stiff-  Figure  5 .  Flow  chart  of  the  wire  icing  program. showing 
ness,  time  stepduration,  and  total  simulation  the  main  program,  where  the  variables  are  initialized 
time.  If  a  droplet  size  distribution  is  spcci-  and  the  results  are  output,  and  the  subroutines. 
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Figure  6  Finite-element  mesh  for 
the  n  ire  icing  program. 


fied,  subroutine  DRPDST  is  called  to  determine  specific- 
droplet  size  information.  Subroutine  FINITE  calculates  the 
air  velocity  field  in  the  vicinity  of  the  wire.  The  finite-element 
grid  for  this  mode)  was  developed  using  MESH  (Albert  and 
Warren  1987),  which  is  an  automatic  finite-element  mesh 
generator.  The  mesh  contains  1 70  nodes  and  288  elements  and 
has  a  half-bandwidth  of  17  (Fig.  6).  The  nodes  are  numbered 
circumferentially  starting  from  the  outside  boundary  of  the 
mesh  and  ending  with  the  nodes  on  the  wire-ice  surface.  Only 
the  windward  half  of  the  wire  is  modeled,  and  the  top,  bottom 
and  upstream  boundaries  of  the  mesh  are  taken  far  enough 
from  the  wire  ( 10  wire  radii)  that  the  flow  velocity  there  is  not 
affected  by  the  wire  and  takes  on  its  uniform  far-fieid  value. 
Two  data  files  are  necessary  to  describe  the  grid.  The  first  file 
contains  the  .v  and  y  coordinates  of  each  node,  and  the  second 
file  is  the  incidence  list,  which  relates  the  element  number  to 
the  three  nodes  in  counterclockwise  order  that  define  the 
element. 

Laplace’s  equation  (eq  1 )  describes  the  air  flow  around  the 
wire.  The  method  of  weighted  residuals  and  Galerkin’s  for¬ 
mulation  were  used  to  develop  the  equations  comprising  the 
local  matrix  for  each  element.  The  subroutine  calculates  the 
terms  of  the  local  matrix  for  each  element  and  assembles  the 
global  matrix  from  the  individual  terms  of  the  local  matrices 
as  they  are  being  generated.  The  global  matrix  is  relatively 


a.  Undeformed  mesh.  h.  Deformed  mesh  after  ice  has  accreted. 

Figure  7.  Finite-element  mesh  dose  to  the  wire. 
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sparse  and  is  stored  in  banded  form.  The  boundary  con¬ 
ditions  are  then  enforced,  and  subroutine  SOLVER  is 
called  to  solve  the  system  of  equations  represented  by 
the  global  matrix.  SOLVER  is  acommercially  available  >. 
routine  that  solves  a  system  of  equations  in  asymmetric  ,§ 
banded  form  using  a  decomposition  procedure .  Finally  t 
the  .vandy  components  of  the  air  velocity  are  determined  g 

for  each  element  using  eq  2.  t> 

Subroutine  FINITE  incorporates  a  movable  bound- 
ary  adjacent  to  the  wire  to  allow  for  the  accretion  of  ice 
on  the  wire.  Figure  8  shows  how  the  boundary  nodes  and 
the  first  three  rows  of  the  finite-element  grid  move  as  the 
ice  grows.  The  mesh  can  deform  to  accept  an  ice  accre¬ 
tion  thickness  up  to  about  1.5  times  the  wire  radius. 

Since  ice  growth  causes  a  change  in  the  air  flow  around  Fi^re  8  Companson  of  lhe  w  col. 
the  accretion,  a  new  velocity  field  is  calculated  during  ,eaion  cfflaency  calculated  by  the 
each  time  step.  present  model  w  ith  that  ofLozowski  and 

Subroutine  TRAJCT  determines  the  point  of  impact  Oleskiw  ( 1983 ). 
on  the  w  ire  of  water  droplets  in  the  air  stream.  This  is 

accomplished  by  starting  with  a  droplet  at  the  right  boundary  of  the  finite-element  grid,  applying  the 
equation  of  droplet  motion  (eq  4),  using  the  drag  coefficient  calculated  from  eq  5.  and  calculating  a 
new  position  and  velocity  of  the  droplet  0. 1  dimensionless  time  units  later.  By  repeating  this  procedure, 
the  droplet  is  tracked  until  it  collides  w  ith  the  wire  or  completely  misses  it.  A  droplet  is  then  inserted 
at  a  new  starting  location  0.025  units  (1  unit=  1  cylinder  radius)  below  the  starting  location  of  the  first 
droplet,  and  this  droplet  is  tracked  until  it  collides  with  the  wire.  This  procedure  is  repeated  while 
moving  the  starting  location  down  the  upstream  boundary  until  the  droplet  misses  the  bottom  of  the 
w  ire. 

Subroutine  TRAJCT  periodicallycallson  subroutine  VELCTY  to  determine  the  airstream  velocity 
at  the  current  droplet  position.  Subroutine  VELCTY  determines  which  finite  element  contains  the 
droplet  by  using  the  following  optimized  technique.  The  subroutine  first  determines  whether  the 
droplet  is  in  the  same  element  as  it  was  previously.  If  not,  it  checks  the  other  elements  by  testing  first 
w  hether  the  droplet  is  within  an.v  andy  tolerance  of  it.  If  it  is,  it  calculates  the  areas  of  three  subtriangles 
formed  by  the  droplet  position  and  two  nodes  of  the  element  in  question.  If  each  of  these  areas  is 
positive,  the  droplet  is  contained  in  that  finite  element.  Using  the  velocities  calculated  in  subroutine 
FINITE,  the  ,t  and  y  components  of  the  air  velocity  at  the  new  droplet  position  are  determined. 

TRAJCT  produces  a  set  of  up  to  80  pairs  of  numbers,  each  pair  consisting  of  Y,  the  ordinate  of  the 
droplet  starting  trajectory,  and  L,  the  distance  along  the  wire  surface  from  the  impact  point  to  the 
stagnation  point  (Fig.  2).  Once  the  path  of  a  particular  droplet  is  calculated,  the  subroutine  finds  the 
point  of  impact  of  the  droplet  on  the  wire-ice  surface.  The  first  step  is  to  determine  whether  the  droplet 
is  past  the  wire  stagnation  point,  and  therefore  within  range  of  a  potential  impact,  or  above  the  top  or 
below  the  bottom  of  the  wire,  and  therefore  cannot  collide  with  the  wire.  The  subroutine  then 
determines  which  segment  of  the  wire-ice  surface  the  droplet  is  near  by  comparing  the  .v  and  y 
coordinates  of  the  wire-ice  boundary  nodes  with  the  coordinates  of  the  current  droplet  position.  If  the 
droplet  has  collided  with  the  wire,  the  intersection  of  the  wire-ice  surface  segment  and  the  last 
trajectory  line  segment  is  determined.  If  not,  another  iteration  is  performed  to  move  the  droplet  closer 
to  the  wire.  Finally ,  L  is  calculated  by  summing  the  segment  lengths  between  the  impact  point  and  the 
stagnation  point.  This  procedure  is  repeated  until  a  set  of  points  is  generated  relating  Y.  the  starting 
ordinate,  and  L.  the  distance  along  the  wire-ice  surface  from  the  stagnation  point  to  the  impact  point, 
for  up  to  80  starting  points.  This  information  is  passed  on  to  subroutine  LSPOLY. 

LSPOLY  is  a  commercial  program  that  uses  a  least-squares  routine  to  fit  the  set  of  data  points 
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obtained  from  TRAJCT  to  a  third-order  polynomial.  The  coefficients  of  the  cubic  equation  are  passed 
on  to  subroutine  COLEFF.  Subroutine  COLEFF  uses  the  cubic  equations  describing  the  trajectory  data 
to  calculate  the  local  collection  efficiency  along  the  wire-ice  surface.  The  starting  ordinate  associated 
with  each  of  the  boundary'  nodes  is  determined  from  the  cubic  relationship,  and  eq  7  in  difference  form 
is  used  to  calculate  the  collection  efficiency  for  each  surface  segment.  Figure  8  is  a  comparison  of  the 
local  collection  efficiency  calculated  in  the  present  model  with  that  of  Lozowski  and  Oleskiw  ( 1 983). 
The  agreement  is  excellent  except  near  the  stagnation  point.  The  difference  is  probably  due  to  the 
discrete  numerical  solution  for  the  air  flow  field  in  this  model.  The  local  collection  efficiencies,  along 
with  the  wind  speed  and  the  liquid  water  content,  are  used  to  calculate  the  mass  flux  of  water  to  each 
surface  segment.  The  mass  of  ice  accreted  on  each  segment  is  determined  by  considering  the 
thermodynamic  processes  occurring  at  the  surface  of  the  accretion  in  subroutine  THERMO. 

THERMO  determines  how  much  of  the  impinging  water  mass  freezes.  The  computer  code  used  in 
THERMO  was  adapted  from  the  model  developed  by  Ackley  and  Templeton  ( 1 979),  which  relied  on 
a  heat  balance  at  the  freezing  surface  to  determine  whether  all,  a  fraction,  or  none  of  the  impinging 
water  freezes.  The  primary  difference  between  the  thermodynamic  model  developed  by  Ackley  and 
Templeton  (1979)  and  this  model  is  that  the  former  treats  the  heat  balance  globally,  while  here  it  is 
calculated  locally.  THERMO  analyzes  the  heat  balance  at  each  of  the  16  surface  segments,  calling 
subroutine  HEATRN  todetermine  the  convective  heat  transfercoefficient.  The  amount  of  ice  accreted 
also  depends  on  the  mass  flux  of  water  to  the  surface  segment,  the  time  step  duration,  the  surface 
segment  length,  and  the  fraction  of  the  incident  water  mass  that  freezes  and  is  calculated  using  eq  1 1 
and  1 2.  Higher  collection  efficiencies  imply  more  incident  water  and  a  greater  heat  flux  to  the  wire- 

ice  surface.  This  results  in  a  surface  temperature  pattern  of 
wanner  temperatures  near  the  stagnation  point  than  at  the  top  and 
bottom  of  the  wire  (Fig.  9).  This  variation  in  w  ire-ice  surface 
temperature  affects  the  density  and  thus  the  volume  of  the  ac¬ 
cretion.  Also,  by  calculating  the  surface  temperature  distribution, 
the  transition  from  dry  grow  th  to  wet  growth  can  be  investigated. 
Since  the  maximum  local  collection  efficiency  of  the  w  ire-ice 
surface  occurs  near  the  stagnation  point,  one  would  expect  the 
transition  to  wet  growth  tooccur  there  first.  By  performingamass 
balance  of  the  runback  water  that  is  created  during  wet  growth 
conditions,  this  phenomenon  could  be  modeled.  Lozowski  et  al. 
(1983)  included  freezing  of  runback  water  in  the  icing  of  a  non¬ 
rotating  cylinder.  Their  model  assumes  symmetric  ice  growth 
about  the  stagnation  line,  which  is  appropriate  for  the  icing  of  a 
non-rotating  cylinder  subject  to  high  winds.  For  the  lower  wind 
speeds  and  finite-stiffness  wires  considered  in  this  model,  deter¬ 
mining  the  direction  of  flow  of  the  unfrozen  water  is  more 
complicated  and  remains  for  future  work. 

The  profile  of  the  accreted  layer  is  calculated  in  subroutine 
PROFIL,  which  calls  on  subroutine  ICEDEN  to  determine  the 
density  of  the  ice  (using  eq  16  and  1 7)  in  each  of  the  sectors  of 
the  accretion  layer.  Subroutine  PROFIL  determines  the  thick¬ 
ness  of  ice  accreted  on  each  surface  segment  during  a  time  step. 
The  accretion  thickness  varies  along  the  wire-ice  surface  because 
of  variations  in  local  collection  efficiency,  fraction  of  incident  water  frozen,  and  ice  density.  It  was 
assumed  that  growth  is  perpendicular  to  the  orientation  of  each  of  these  segments,  so  subroutine 
PROFIL  distributes  the  mass  accreted  on  each  surface  segment  in  a  trapezoidal  cross-sectional  shape, 
as  shown  in  Figure  10a.  The  coordinates  of  these  trapezoids  are  then  passed  to  subroutine  ROTATE, 
where  the  torque  due  to  the  asymmetric  ice  weight  is  calculated. 
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Figure  9.  Example  of  wire -ice 
surface  temperatures  calculated 
by  subroutine  THERMO  for  an 
air  temperature  of  ~15°C. 
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a.  Ice  layer  accreted 
in  each  sector. 


c.  New  ice  profile  in 
each  sector. 


Figure  10.  Steps  in  the  c 
ice  layer  profile. 


h.  Wire  rotated  due  to 
eccentric  ice  mass. 


d  Smoothed  ice  ac¬ 
cretion  profile. 


'ale  ulation  of  the  accreted 


The  purpose  of  subroutine  ROTATE  is  to  calcu¬ 
late  the  torque  created  by  the  accretion  of  ice  on 
the  windward  side  of  the  wire,  the  amount  of  ro¬ 
tation  due  to  this  torque,  and  the  coordinates  of 
the  finite-element  boundary  nodes  after  rotation. 
Using  eq  20,  ROTATE  calculates  the  incremen¬ 
tal  rotation  at  the  end  of  each  time  step  (Fig.  10b). 
The  post-rotation  cross-sectional  area  of  the  new 
ice  accretion  layer  in  each  sector  is  computed  and 
distributed  in  a  trapezoidal  shape,  as  was  done  in 
PROFIL  (Fig.  10c).  ROTATE  then  calculates  a 
continuous  accretion  boundary  by  averaging  the 
adjacent  end  point  coordinates  of  each  pair  of  the 
new  sector  boundaries  (Fig.  lOd). 

RESULTS 

The  results  of  13  icing  simulations  are  de¬ 
scribed  in  this  section.  The  simulations  investi¬ 
gate  the  variation  in  vertical  (gravitational)  and 
horizontal  (wind)  load  on  a  3.5-cm-diameter 
wire  due  to  variations  in  the  atmospheric  con¬ 
ditions  and  the  transmission  line  design.  The 
simulations  are  carried  out  for  one  hour.  A  typ¬ 
ical  icing  event  is  represented  by  the  default  case 
for  which  U0  =  10  m/s,  7>  -1 5°C.  IV  =  0.0005 
kg/m-’  and  r  =  1 7  pm.  The  default  wire  properties 
are  taken  as  D  =  3.5  cm.  Lw=250  m  and  K  =  1 6 


Table  1.  Results  of  icing  simulations.  For  each  simulation  the  atmospheric  and  wire  variables  are 
given  as  well  as  the  final  wire  rotation  angle,  the  accreted  mass  per  unit  length  of  wire,  the  vertical  load 
on  the  wire  span,  the  horizontal  load  on  the  wire  span,  and  the  accretion  diameter  perpendicular  to 
the  wind.  The  default  case  isshown  first.  In  the  other  cases  the  parameter  values  that  differ  from  those 
in  the  default  case  are  in  bold  type. 
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Nm/rad.  The  final  cross-sectional  accretion  diameter 
perpendicular  to  the  wind,  the  rotation  angle,  the  mass 
of  accretion  per  meter  of  wire,  the  ice  load  and  the 
wind  load  are  show  n  in  Table  1  for  all  13  simulations. 
The  w ind  load  is  calculated  using  eq  20  with  Cd  =  1 . 
The  accretion  profiles  are  compared  in  Figure  1 1. 

Atmospheric  variables 

The  ambient  air  temperature  was  varied  from  the 
default  value  of  -1 5°C  to  -5  and  -25°C.  The  effect  of 
these  different  air  temperatures  is  to  change  the  den¬ 
sity  of  the  accretion,  which  alters  the  accretion  profile 
and  thus  the  moment  arm  causing  rotation.  At  the 
higher  air  temperatures,  the  surface  temperature  of  the 
accretion  approaches  0°C  near  the  stagnation  point 
and  not  all  the  incoming  liquid  water  freezes.  Assum¬ 
ing  this  water  is  shed  rather  than  running  back  and 
freezing  gives  a  lower  accretion  mass  than  in  the 
colder  simulations.  Table  1  indicates  that  higher  am¬ 
bient  temperatures  lead  to  slightly  lower  ice  and  wind 
loads  on  the  transmission  line.  The  accretion  profiles 
are  shown  in  Figure  1  la. 

The  wind  speed  was  varied  from  the  default  value 
of  1 0  m/s  to  20  and  30  m/s.  The  increased  wind  speed 
increases  the  collection  efficiency  of  the  wire  and  the 
flux  of  water  droplets  past  the  wire,  resulting  in 
greatly  increased  ice  and  wind  loads  on  the  w  ire.  The 
low-wind-speed  accretion  has  a  lower  surface  tem¬ 
perature  and  is  likely  to  be  rougher  than  the  ice  accret¬ 
ed  at  the  higher  wind  speeds.  The  accretion  roughness 
will  also  affect  the  wind  load  on  the  wire,  but  this  effect 
has  not  been  quantified.  The  accretion  profiles  are 
show  n  Figure  1 1  b. 

The  liquid  water  content  was  varied  from  its  de¬ 
fault  value  of  0.0005  kg/m ''  to  0.00025  and  0.00075 
kg/nv\  The  ice  load  is  roughly  proportional  to  the 
liquid  water  content  for  these  simulations.  The  w  ind 
load  increases  somewhat  with  liquid  water  content 
because  of  the  associated  increase  in  the  cross-sec¬ 
tional  diameter  of  the  accretion.  The  accretion  profiles 
are  shown  in  Figure  1  lc. 

The  droplet  radius  was  varied  from  its  default 
value  of  1 7  pm  to  1 0  and  1 3  pm.  The  decreased  col¬ 
lection  efficiency  for  the  smaller  droplet  sizes  gives 
smaller  ice  loads  and  somewhat  smaller  wind  loads. 
The  accretion  profiles  are  shown  in  Figure  1  Id. 

Wire  variables 

The  average  wire  stiffness  was  varied  from  its  de¬ 
fault  value  of  1 6  Nm/rad  to  8  and  32  Nm/rad.  This  was 


a.  Air  temperatures  -5.  -15  and 
-25  T. 


X 


h.  Wind  speeds  10.  20  and  30  mis. 


\ 


. -  0  75  g  m3 

e.  Liquid  water  contents  0.25.  0.5 
and  0.75  glnr 

Figure  1 1 .  Final  ice  accretion  profiles  for  15 
model  simulations.  For  the  default  case,  air 
temperature  =  -l5°C.  wind  speed  =  10  mis. 
liquid  water  content  =  0.5  glnr1,  drop  size  = 
17  pm,  wire  stiffness  =  16  Nm/rad,  wire 
length  =  250  m.  For  the  oilier  cases  one 
variable  at  a  time  was  changed.  Each  figure 
shows  the  default  case  along  with  two  other 
cases. 
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assumed  to  have  been  done  by  using  different 
materials  in  manufacturing  the  wire,  without 
changing  the  wire  diameter.  The  wire  stiffness 
affects  the  amount  of  rotation  for  a  given  accreted 
mass  of  ice.  which  affects  the  accretion  collection 
efficiencies  and  the  final  accretion  profile.  The  ice 
and  wind  loads  decrease  slightly  as  the  wire  stiff¬ 
ness  increases.  The  accretion  profiles  are  shown  in 
Figure  1  le. 

Changing  the  span  of  the  transmission  line 
between  support  structures  changes  the  average 
stiffness  as  well  as  the  length  of  the  wire  (eq  29). 
The  wire  length  was  varied  from  its  default  value 
of  250  m  to  125  and  500  m.  Assuming  the  same 
w ire  shear  modulus  and  polar  moment  of  inertia 
gives  wire  stiffnesses  of  8  and  32  Nm/rad.  re¬ 
spectively.  for  these  cases.  The  longer,  more- 
flexible  wire  accretes  more  ice  and  rotates  through 
a  larger  angle  than  the  stiffer  wires.  This  effect  is 
compounded  in  the  ice  and  w  ind  load  calculations 
by  the  greater  w  ire  length  over  which  the  ice 
collects,  resulting  in  greatly  increased  vertical  and 
horizontal  loads  on  support  structures.  The  accre¬ 
tion  profiles  are  shown  in  Figure  1  If. 

Summary 

These  simulations  quantify  the  dependence  of 
ice  and  wind  loads  on  atmospheric  varibles  and  on 
design  \  ariables  of  the  wire  and  support  structures. 
Ice  and  wind  loads  on  a  transmission  line  increase 
w  ith  increasing  w  ind  speed,  liquid  water  content 
and  drop  size.  The  wind  load  increases  somewhat 
with  decreasing  airtcmperature.  but  the  ice  load  is 
not  sensitive  to  air  temperature  variations.  Ice  and 
wind  loads  decrease  with  increasing  wire  stiffness 
but  only  slightly.  However,  these  loads  increase 
markedly  with  increasing  wire  length  and  the 
associated  decrease  in  average  wire  stiffness. 

CONCLUSIONS 
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d.  Drop  sizes  10,13  and  17  pm. 


e.  Wire  stiffnesses  8.  16  and  32  N mi  rad. 


f.  Wire  lengths  125.  250  and  500  m  with 
corresponding  stiffnesses. 

Figure  11  (cont'd).  Final  ice  accretion  profiles 
for  13  model  simulations.  For  the  default  case, 
air  temperature  =  -15X1,  wind  speed  =  10  m/s, 
liquid  water  content  =  0.5  g/nr,  drop  size  =  17 
pm,  wire  stiffness  =16  Nm/rad,  wire  length  = 
250 m.  For  the  other  cases  one  variable  at  a  rime 


After  identifying  a  potential  or  existing  ice  ,V<M  changed.  Each  figure  shows  the  default 
accretionproblem.thedesignengineermustchoose  case  a^on8  'v'th  two  other  tases. 
between  designing  for  a  particular  icing  event  at 
the  location  of  interest  or  relocating  the  transmission  line. 

Transmission  line  relocation  should  always  be  the  first  option  investigated  because  avoiding  a 
foreseeable  problem  is  the  easiest  way  to  solve  it.  This  may  be  an  attractive  option  if  the  line  traverses 
terrain  where  atmospheric  icing  is  a  problem  for  only  a  short  distance,  such  as  a  mountain  pass.  It  may 
be  possible  to  adjust  the  transmission  line  orientation  and  location  in  a  pass  to  avoid  the  full  impact 
of  prevailing  winds.  The  wind  speed  has  a  great  effect  on  the  ice  load  and  the  wind  load  because  of 


the  increase  in  collection  efficiency  at  higher  wind  speeds,  the  increased  flux  of  droplets  past  the  wire, 
and  the  dependence  of  wind  drag  on  the  square  of  the  wind  speed.  In  many  cases,  however,  right-of- 
way  acquisition  may  preclude  relocation  of  the  line.  The  transmission  line  may  also  traverse  open  and 
exposed  areas  for  miles,  in  which  case  relocation  is  simply  not  an  option.  If  the  problem  cannot  be 
avoided,  the  design  of  the  line  and  the  support  structures  must  either  take  into  account  the  increased 
loading  associated  with  atmospheric  icing  or  bear  the  losses  associated  with  increased  maintenance 
associated  with  icing  storms  and  the  loss  of  revenue  from  customers  during  power  outages.  If  a  severe 
icing  storm  occurs  only  once  in  20  years,  this  may  be  an  attractive  option. 

The  computer  simulations  indicate  that,  with  regard  to  structural  design,  the  length  or  stiffness  of 
a  transmission  line  can  be  changed  to  decrease  the  ice  and  wind  loads  on  the  wire  and  support  structures 
due  to  an  icing  storm.  The  wire  stiffness  can  be  increased  in  practice  by  increasing  the  diameter  of  the 
wire  (which  also  decreases  the  wire  collection  efficiency),  by  increasing  the  shear  modulus  of  the 
material  in  the  wire,  or  by  using  phase  spacers.  Phase  spacers  are  short  rigid  cylinders  made  from  an 
insulating  material  such  as  fiberglass  and  are  installed  midspan  with  a  primary  purpose  of  preventing 
two  phases  of  the  transmission  system  from  touching  and  shorting  out.  These  spacers  can  be  designed 
to  prevent  the  wire  from  rotating  where  they  attach,  which  effectively  increases  the  average  wire 
stiffness.  An  additional  advantage  is  that  these  spacers  can  be  installed  on  transmission  lines  already 
in  service.  It  must  be  emphasized  that  model  results  indicate  only  a  5-10%  reduction  in  the  ice  load 
or  w  ind  load  on  the  wire  and  support  structures  associated  with  doubling  the  average  w  ire  stiffness. 
However,  the  reduction  in  support  structure  ice  and  wind  loads  associated  with  reducing  the  wire 
length  (which  also  causes  the  average  wire  stiffness  to  increase)  is  roughly  the  same  as  the  length 
reduction.  This  effect  is  primarily  due  to  the  decreased  length  of  wire  between  supports.  In  the  wire 
simulations  presented  earlier,  for  example,  a  reduction  in  the  wire  length  from  250  m  to  1 25  m  caused 
a  53%  reduction  in  both  the  ice  and  wind  loads  on  the  support  structures. 

These  results  are  from  a  limited  numba  of  simulations  for  which  the  duration  of  the  icing  event  was 
only  one  hour,  while  many  icing  storms  last  hours  or  even  days.  They  require  confirmation  through 
field  or  wind  tunnel  testing  before  confident  design  changes  can  be  made  for  planned  or  existing 
transmission  lines.  It  must  also  be  considered  that  changing  the  w  ire  stiffness  or  support  structure 
separation  may  not  be  practical  from  the  standpoint  of  designing  the  line  for  its  primary  purpose. 

FUTURE  WORK 

Future  work  can  be  divided  into  three  areas:  modeling,  laboratory  studies  and  field  studies.  The 
limitations  of  modeling  are  the  physical  simplifications  that  must  be  made  to  quantify  many  of  the 
factors  that  contribute  to  the  complex  phenomenon  of  atmospheric  ice  accretion.  Laboratory  studies 
are  limited  by  the  sophisticated  hardware  that  must  be  employed  to  simulate  an  icing  event.  Field 
studies  are  difficult  because  of  the  harsh  conditions  associated  with  icing  storms  and  the  infrequency 
of  their  occurrence.  However,  more  laboratory  and  field  work  is  required  to  verify  results  obtained  in 
this  and  other  similar  computer  modeling  efforts. 

Future  modeling  studies  can  be  improved  in  two  ways:  improving  the  thermodynamic  calculations 
at  the  surface  of  the  accretion  and  improving  the  calculation  of  the  wind  load.  Many  of  the  ideas 
presented  in  this  study  require  verification,  especially  the  concept  of  the  accretion  having  a  variable 
temperature  distribution  on  its  surface.  A  method  of  quantifying  the  change  in  convective  heat  transfer 
coefficient  along  the  surface  of  an  accretion  should  be  investigated.  To  improve  the  accuracy  of  wind 
load  calculations,  future  work  should  include  a  parameterization  of  the  accretion  roughness  and  its 
effect  on  the  drag  coefficient.  Future  work  might  also  include  an  analysis  of  the  fate  of  the  unfrozen 
water  on  the  surface  of  a  glaze  ice  accretion. 

Laboratory  studies  would  be  useful  in  verifying  the  results  presented  in  this  study,  especially  the 
variation  in  temperature  on  the  surface  of  the  accretion.  Wind  tunnel  experiments  could  be  set  up  to 
measure  the  variation  in  convective  heat  transfer  coefficient  on  the  surface  of  a  cylinder  with  a 
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noncylindrical  cross  section  to  extend  the  work  of  Achenbach  (1974).  Measurements  of  the  drag 
coefficient  associated  with  different  roughness  heights  would  also  be  useful. 

Finally,  assumptions  and  conclusions  should  be  verified  in  the  field.  Continued  measurements  of 
ice  accretion  and  w  ind  loads  on  actual  transmission  lines  by  Ontario  Hydro  Research  (Krishnasamy 
1983)  would  be  of  great  value.  In  the  final  analysis,  field  inspections  and  measurements  by  technical 
and  nontechnical  personnel  will  probably  contribute  the  most  to  our  understanding  of  the  ice  accretion 
process  and  improved  engineering  design. 
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APPENDIX  A:  DETERMINING  THE  ANGLE  OF  TWIST  ALONG  A  CYLINDER 
FIXED  AT  BOTH  ENDS  SUBJECT  TO 
A  DISTRIBUTED  TORQUE  ALONG  ITS  LENGTH 


First  determine  the  angle  of  twist  along  a  cylinder  fixed  at  one  end  and  subject  to  a  concentrat¬ 
ed  torque  at  r  =  L: 


'A 


GJ 


From  Volterra  and  Gaines  (1971,  p.  217)  the  angle  of  twist  at  a  distance  z  from  the  fixed  end  is 

0(c)  =  :T/GJ  (Al) 

where  G  =  shear  modulus 

./  =  polar  moment  of  inertia 
T  =  applied  torque. 

The  torsional  stiffness  K  of  the  cylinder  is  the  ratio  between  the  applied  torque  and  the  rotation 
angle 


K  =  GJ/:. 


(A2) 


Next  determine  the  angle  of  twist  along  a  cylinder  fixed  at  both  ends  subject  to  a  concentrated 
torque  at  :  =  I. 


T 


The  resisting  moments  at  the  fixed  ends  must  satisfy  7j+7\  =  T.  At  e  =  /  the  angles  of  twist  must 
match:  0|  =  7j//G7  must  equal  02  =  7\(Z.-/)/GJ. 

Solving  for  7j  gives 

7j  =  {L-DTIL  (A3) 


and 


0(/>  =  (L-l)ITIGJl  (A4) 

How  does  the  twist  angle  depend  on  r?  It  must  vary  linearly  on  either  side  of  the  concentrated 
torque  as  in  eq  Al  and  be  zero  at  the  fixed  ends,  so 
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0(r)  =  (L-l)zT/GJL 

z<l 

0(r)  =  l(L-z)T/GJL 

z>  !. 

(A5) 


Now  determine  the  angle  of  twist  along  a  cylinder  fixed  at  both  ends  subject  to  a  distributed 
torque  t  per  unit  length  along  the  whole  length.  First  look  at  a  small  element  dr: 


The  angle  of  twist  at  r  due  to  the  applied  torque  /  over  the  length  db  is  (using  eq  A4): 

d0(r)  =  (L-z)ztdz/GJL.  (A6) 

The  twist  angle  at  any  other  point  z0  due  to  that  applied  torque  is  (using  eq  A5) 


d0(-o)  =  (L-z)z0rdz/GJL 
d0(r„)  =  :{L-:0)tdz/GJL 


(A7) 


The  distributed  torque  I  per  unit  length  is  applied  over  the  whole  length  of  the  cylinder,  so  eq  A7 
can  he  integrated  over  r  to  get  the  total  twist  angle  at  z0: 


0(_-  )  =  [' "" iL-  :Jil:  dr  +  f  it 

G.IL  OJL 


(AS) 


vv  Inch  gives 

0(:„>  =  vL)U\-:/L )  ( zIDIlG.I .  (A9) 

This  is  eq  21  in  the  main  text. 
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